#include <iostream>
#include <iomanip>
#include <array>
#include "Core/numlib.h"
#include "Core/TensorExpr.h"
#include "Embedded/EmbeddedGrid.h"
#include "Embedded/EmFuncFiller.h"
#include "Embedded/EmLevelOp.h"
#include "Embedded/PartialDer.h"
#include "../example/common/Ellipsoid.h"
#include "../example/common/TestUtility.h"

template class Tensor<Real, SpaceDim>;

const Real gTol = 1e-12;

template <int Dim, int pD, int desiredOrder, class T1, class T2>
std::array<Real,3> testPartialDer(int N, const T1 &exactSol, const T2 &exactGrad)
{
  using rVec = Vec<Real, Dim>;
  const int numGhost = desiredOrder / 2;
  RectDomain<Dim> rd(Box<Dim>(0, N-1), 1.0/N, CellCentered, numGhost);
  Ellipsoid<Dim> ellipsoid({0.5, 0.5, 0.5}, {0.25, 0.125, 0.25});
  EmbeddedGrid<Dim> eGrid(rd, [&](rVec q, rVec &closest, rVec &normal, bool &contained, Real gTol) {
    ellipsoid.check(q, closest, normal, contained, gTol);
    normal = -normal;
    contained = !contained;
  }, true, 0.2, gTol);
//#ifndef NDEBUG
//  if(get_dbglevel() >= 3)
//    eGrid.dump(tmpos);
//#endif
  PartialDer<Dim, pD, desiredOrder+1> pd(eGrid);

  EmFuncFiller<Dim> eff(eGrid);
  EmLevelOp<Dim>    elv(eGrid);
  Box<Dim> gbx = rd.getGhostedBox();
  Tensor<Real, Dim> phi(gbx);
  Tensor<Real, 1> bcData(eGrid.numGhosts);
  eff.fillPointValues(phi, CellCentered, exactSol);
//  eff.fillBdryValues(phi, "curvilinear", exactSol);
  eff.fillBdryValues(bcData, "curvilinear", exactSol);
  eff.transferBdryValues(phi, bcData);
  Tensor<Real, Dim> err(gbx);
  eff.fillPointValues(err, CellCentered, exactGrad);
  Tensor<Real, Dim> numpd(gbx);
  pd.apply(phi, numpd);
  elv.zeroExtNodes(numpd);
  err = numpd - err;
  elv.zeroExtNodes(err);

  std::array<Real, 3> errnorm;
  for(int p = 0; p <= 2; ++p)
    errnorm[p] = elv.NRM(err, CellCentered, p);
  return errnorm;
}

template <int Dim, int pD, int desiredOrder, class T1, class T2>
void doTest(const std::vector<int> &gridSize, const T1 &exactSol, const T2 &exactGrad)
{
  const int numGrid = gridSize.size();
  std::vector<std::array<Real, 3>> errnorm;
  for(int n=0; n<numGrid; ++n) {
    std::cout << "Testing h = " << oneOver(gridSize[n]) << std::endl;
    errnorm.push_back(testPartialDer<Dim, pD, desiredOrder>(gridSize[n], exactSol, exactGrad));
  }
  printConvergenceTable(&gridSize[0], errnorm);
}

int main()
{
  reset_dbglevel(0);
  const int pD = 0;
  std::cout << "SpaceDim = " << SpaceDim << std::endl;
  std::cout << "========================================================" << std::endl;
#if DIM == 2
  auto exactSol = [](const Vec<Real, 2> &x) {
    return sin(2 * M_PI * x[0]) * cos(M_PI * x[1]);
  };
  auto exactGrad = [](const Vec<Real, 2> &x) {
    return 2 * M_PI * cos(2 * M_PI * x[0]) * cos(M_PI * x[1]); // pD = 0
    //return -M_PI * sin(2*M_PI*x[0]) * sin(M_PI*x[1]); // pD = 1
  };
  std::cout << "Testing fourth-order ..." << std::endl;
  doTest<SpaceDim, pD, 4>({256, 128, 64, 32}, exactSol, exactGrad);
  std::cout << "Testing sixth-order ..." << std::endl;
  doTest<SpaceDim, pD, 6>({256, 128, 64, 32}, exactSol, exactGrad);
  std::cout << "Testing eighth-order ..." << std::endl;
  doTest<SpaceDim, pD, 8>({128, 64, 32}, exactSol, exactGrad);
#elif DIM == 3
  auto exactSol = [](const Vec<Real, 3> &x) {
    return sin(M_PI * x[0]) * cos(M_PI * x[1]) * sin(M_PI * x[2]);
  };
  auto exactGrad = [](const Vec<Real, 3> &x) {
    return M_PI * cos(M_PI * x[0]) * cos(M_PI * x[1]) * sin(M_PI * x[2]); // pD = 0
  };
  std::cout << "Testing fourth-order ..." << std::endl;
  doTest<SpaceDim, pD, 4>({128, 64, 32, 16}, exactSol, exactGrad);
  std::cout << "Testing sixth-order ..." << std::endl;
  doTest<SpaceDim, pD, 6>({128, 64, 32, 16}, exactSol, exactGrad);
  std::cout << "Testing eighth-order ..." << std::endl;
  doTest<SpaceDim, pD, 8>({128, 64, 32, 16}, exactSol, exactGrad);
#endif
  return 0;
}

// The linear system for (20,29,43) has condition number 9.21e+04.
